home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / invhess_r < prev    next >
Encoding:
Text File  |  1994-12-20  |  1.8 KB  |  63 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Inverse of an upper Hessenberg matrix.
  4.  
  5. // Syntax:      invhess ( X , Y )
  6.  
  7. // Description:
  8.  
  9. //      invhess(x, y), where X is an N-vector and Y an N-1 vector, is
  10. //      the matrix whose lower triangle agrees with that of
  11. //      ones(N,1)*X' and whose strict upper triangle agrees with that
  12. //      of [1, Y]*ONES(1,N). 
  13.  
  14. //      The matrix is nonsingular if X[1] != 0 and X[i+1] != Y[i] for
  15. //      all i, and its inverse is an upper Hessenberg matrix. If Y is
  16. //      omitted it defaults to -X[1:N-1]. 
  17.  
  18. //      Special case: if X is a scalar invhess(X) is the same as
  19. //      invhess(1:X). 
  20.  
  21. //      References:
  22. //       F.N. Valvi and V.S. Geroyannis, Analytic inverses and
  23. //            determinants for a class of matrices, IMA Journal of Numerical
  24. //            Analysis, 7 (1987), pp. 123-128.
  25. //       W.-L. Cao and W.J. Stewart, A note on inverses of Hessenberg-like
  26. //            matrices, Linear Algebra and Appl., 76 (1986), pp. 233-240.
  27. //       Y. Ikebe, On inverses of Hessenberg matrices, Linear Algebra and
  28. //            Appl., 24 (1979), pp. 93-97.
  29. //       P. Rozsa, On the inverse of band matrices, Integral Equations and
  30. //            Operator Theory, 10 (1987), pp. 82-95.
  31.  
  32. //    This file is a translation of invhess.m from version 2.0 of
  33. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  34. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  35.  
  36. //-------------------------------------------------------------------//
  37.  
  38. invhess = function ( x , y )
  39. {
  40.   local (x, y)
  41.  
  42.   n = max(size(x));
  43.   //  Handle scalar x.
  44.   if (n == 1)
  45.   {
  46.     n = x;
  47.     x = 1:n;
  48.   }
  49.  
  50.   x = x[:];
  51.  
  52.   if (!exist (y)) { y = -x; }
  53.   y = y[:];
  54.  
  55.   A = ones(n,1)*x';
  56.   for (j in 2:n)
  57.   {
  58.     A[1:j-1;j] = y[1:j-1];
  59.   }
  60.  
  61.   return A;
  62. };
  63.